
%MACRO propensity (in,out);

*****************************************************************************
* Program: PROPENSITY					   						
* 																			
* Description: Computes propensity scores and associated weights for use in
* 				regression. Performs diagnostic checks (common support,
*				standardized differences and empirical cumulative
*				distribution functions).
*						
* Last modified: September 2, 2015												
* Created by: Martin Wegman													
*****************************************************************************;

%let dat = data; * Dataset name stem;

%let v = 30; * From missing_stage2 sas file;

data loc.&dat.&v;
	set loc.&in;
run;

****************************************************************************;
****************************************************************************;
************************    Initial Processing     *************************;
****************************************************************************;
****************************************************************************;

* Recode site vars to ease analytic interpretation;

proc format library = loc;
value site2_ 1="VTC"
			0="CDDC";
value site3_		1="CDDC"
			0="1VTC";
value site3b_		1="CDDC"
			0="VTC";
value combin_ 	1="SR+SR"
				2="SV14+SR"
				3="SV0+SR"
				4="SR+SV90"
				5="SV14+SV90"
				6="SV0+SV90"
				7="SR+SV30"
				8="SV14+SV30"
				9="SV0+SV30";
run;

%let v2 = %eval(&v + 1);

data loc.&dat&v2;
	set loc.&dat&v;
	if site = 1 then site2 = 0;
	else if site = 2 then site2 = 1;
	
	if site = 1 then site3 = 1;
	else if site = 2 then site3 = 0;

	if site = 1 then do;
		ih05_totaltimes = ih05_totaltimes  -1; * Since this will be used in propensity model, 
													must reflect characteristic prior to group assignment;
	end;
	if site = 2 then do; * Same as above, except add 1 to VTC group to avoid log of 0;
		ih02_log = log(exp(ih02_log) + 1);
	end;

	format site2 site2_. site3 site3_. combination combin_.;
run;

****************************************************************************;
****************************************************************************;
************************     Model variables       *************************;
****************************************************************************;
****************************************************************************;

* Characteristics before treatment is assigned;

* Continuous vars;
%let var1 = 		dm0
					ih01_sqrt
					ih02_log
					ih05_totaltimes 
					du16_log 
					du02b_heroin_years 
					ss_so ss_fam ss_fri; 

* Binary vars;
%let var2b =		dm3_m
					dm4_m
					dm5_m  
					rbd1
					daily_heroin  
					lifetime_alcohol 
					lifetime_other_opiates 
					lifetime_benzos 
					lifetime_stimulants ;

* >2 level categorical vars;
%let var2m = 		dast_cat_m 
					dm2_m;

* Vars after treatment assignment;

* Continuous vars;
%let var3 = 	 cs01_craving
					so_rec_d so_amb_d so_ts_d; 
* Categorical vars;
%let var4 = ;

****************************************************************************;
****************************************************************************;
********************      Compute stabilized IPTWs         *****************;
****************************************************************************;
****************************************************************************;

/* Code adapted from Vaughn et al 2015;

/* First model the exposure of interest (i.e. the treatment) as a function of
potential confounders. */

/* The output dataset ps_p contains the IPTW weights */

%let v = 31;

proc sort data = loc.&dat&v; by combination _imputation_; run;

ods exclude all;
proc logistic data = loc.&dat&v;
	by combination _imputation_;
	class site3 &var2b &var2m;
	model site3 (ref= FIRST) = &var1 &var2b &var2m/
	lackfit outroc = ps_r1a; * ref = FIRST sets CDDC as reference;
	output out= loc.ps_p XBETA=beta STDXBETA= betasd PREDICTED = ps_pred;
run;
ods exclude none;

/* Calculate weights based on predicted probabilities*/
data loc.ps_weight;
	set loc.ps_p;
	if site3 = 1 then wt_ps_pred = ps_pred;
	else wt_ps_pred = 1-ps_pred;
run;

/* Calculate a mean weight by exposure group in order to standardize weights
*/
proc sort data= loc.ps_weight;
	by combination _imputation_ site3;
proc means data = loc.ps_weight noprint;
	var wt_ps_pred;
	by combination _imputation_ site3;
	output out = loc.q mean = mn_wt;
run;

/* Stabilize the weights with the mean for each treatment group*/
/* The variable wt2 is the standardized IPTW weight */
data loc.weights;
	merge loc.q loc.ps_weight;
	by combination _imputation_ site3;
	wt2 = mn_wt/wt_ps_pred;
	drop _type_ _freq_;
run;


data loc.&out;
	set loc.weights;
run;

****************************************************************************;
****************************************************************************;
*********************   Remove old datasets               ******************;
****************************************************************************;
****************************************************************************;


proc datasets lib=loc nolist;       
	delete weights q ps_weight &dat.30 &dat.31 ps_p;   
run;
quit;

****************************************************************************;
****************************************************************************;
*****************                   End                        *************;
****************************************************************************;
****************************************************************************;

%MEND propensity;
